Revisiting spherically symmetric relativistic hydrodynamics 
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In this paper we revise two classical examples of Relativistic Hydrodynamics in order to illustrate 
in detail the numerical methods commonly used in fluid dynamics, specifically those designed to 
deal with shocks, which are based on a finite volume approximation. The two cases we consider are 
the relativistic blast wave problem and the evolution of a Tolman-Oppenheimer-Volkoff star model, 
in spherical symmetry. In the first case we illustrate the implementation of relativistic Euler's 
equations on a fixed background space-time, whereas in the second case we also show how to couple 
the evolution of the fluid to the evolution of the space-time. 

PACS numbers: 95.30.Lz, 02.60.-x 04.20.-q 
I. INTRODUCTION 

In this paper we present the elements needed to im- 
plement the numerical solution to the relativistic Eulcr 
equations in spherical symmetry using examples in spe- 
cial and general relativity. The aim is to describe the 
necessary tools to implement numerical codes able to 
deal with basic problems involving relativistic hydrody- 
namics which eventually are used to model high energy 
astrophysical phenomena, for instance stellar collapse of 
compact stars like supernovae core collapse, shock waves 
propagating out from a spherical source interacting with 
the interstellar medium, etc. We are particularly inter- 
ested in revising the specific numerical methods that are 
commonly used and present a detailed flavor of high reso- 
lution shock capturing methods used in general relativis- 
tic hydrodynamics. We focus on the description of two 
representative physical cases: the spherically symmetric 
blast wave on Minkowski space-time, which is considered 
to be a test case in hydrodynamics and the evolution of 
a Tolman-Oppenhcimcr-Volkoff (TOV) star model, made 
of a self-gravitating polytropic ideal gas on a dynamical 
space-time background. 

We consider the problem of solving relativistic hydro- 
dynamical systems as an initial value problem, ruled by 
relativistic Euler's equations. We have to provide initial 
data for a relativistic fluid, which then evolves according 
to Euler's equations in the blast wave case, and in the 
case of the TOV star we need to solve simultaneously 
Euler's and Einstein's equations. 

In the blast wave case we start up with free initial 
data, which we choose to correspond to an ideal gas dis- 
tributed into two concentric spherical chambers, being 
the inner one where the gas is at high pressure whereas 
in the outer sphere the pressure is smaller. In the TOV 
star case it is not possible to choose arbitrary initial data, 
it is necessary to construct initial data that are consistent 
with Einstein's equations. In order to have a complete 
description of the second problem we introduce the nec- 
essary general relativistic background. 

We base our numerical treatment on an Eulerian de- 



scription of the fluid equations using a flux balance form 
of the system of equations, which requires the definition 
of conservative variables on top of the discrete cell cen- 
tered mesh. The solution of relativistic Euler's equa- 
tions is found using a High Resolution Shock Capturing 
method (HRSC) based on the approximate solution of 
local Riemann problems at the intercell boundaries. Par- 
ticularly we use the HLLE Riemann approximate solver 
with a linear piecewise reconstructor of variables. Being 
this a paper revisiting the numerical methods, we try to 
be as specific as possible in the description of each of the 
steps within the appropriate section. 

The paper is organized as follows. In section |H] we set 
the equations of hydrodynamics for a spherically sym- 
metric space-time, define conservative variables and set 
the system of Euler's equations as a flux balance set of 
equations, so as the numerical methods used for the so- 
lution. In IIIII we present the blast wave case and describe 
in detail its properties and in IIVI we present the evolu- 
tion of TOV stars. Finally in [V] we present some final 
comments. 



II. HYDRODYNAMICS 

A. The equations of relativistic hydrodynamics in 
spherical symmetry 

As a starting point we describe a spherically symmetric 
space-time line element in spherical coordinates to be of 
the form 



ds 2 = -a 2 (t, r)dt 2 + a 2 (t, r)dr 2 + r 2 d9 2 + r 2 sin 2 

(1) 

where t is the time coordinate and (r, 9, cf>) are the usual 
spherical coordinates and where we have assumed geo- 
metric units where G = c = 1. This line-element will 
serve to workout the two problems we deal with: the 
hydrodynamics onto the Minkowski space-time where 
a = a = 1 and the TOV star where such metric func- 
tions obey Einstein's equations. We choose the matter 
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model to correspond to a perfect fluid, which means the 
fluid is not subject to heat transfer and viscosity effects; 
the stress energy tensor of a perfect fluid given in general 
relativistic form reads [l|, [|[ : 



9 t (aar 2 w') H ^d r (aar z v r ) = 



2„.r\ 



(6) 



T^ v = p huV+pg 



(2) 



where g is the determinant of the metric tensor in ([1} and 
therefore y/—g = aar 2 sinO. Defining D = poW = poau* 



where p is the rest mass density of the fluid, p its pres- and considering aa ^ Q we have 
sure, u M is the four velocity of the fluid on the space-time 
described by ([TJ) and h the specific enthalpy 



h = 1 + e + p/po 



(3) 



where e is the specific internal energy of the gas. The 
stress-energy tensor © is commonly found in the litera- 
ture written as T^ v = (p + p)u M u v +pg^, where p is the 
total energy density of the gas which can be expressed 
in terms of the rest mass density and internal energy as 
p = po(l + e). Also g^ v are the components of the inverse 
of the 4- metric in (|T|). 

The fluid equations arc given by the local mass conser- 
vation law and the Bianchi identity, that are respectively: 



V„(pou") = 0, 
V M r" v = 0, 



(4) 
(5) 



where V M is the covariant derivative consistent with ((T|). 
When these equations are projected onto space-like hy- 
persurfaces and their normal directions one obtains the 
relativistic Euler equations [3 5]. The result is a set of 
equations for the primitive variables po,v r ,p or cquiva- 
lently po,v r , e, where v r is the three velocity of the fluid 
elements measured by an Eulerean observer. The way to 
relate the spatial velocity with the spatial components 
of the four velocity of the fluid in <j2j) is using the re- 
lation v r = u r ^l — g rr v r v r = u r Vl — a 2 v r v r = u r /W, 
where W is the Lorentz factor, which in turn is defined 
by W = au l . 

It is well known that Euler's equations develop discon- 
tinuities in the hydrodynamical variables even if smooth 
initial data are considered @ . Therefore one may use as a 
first try a finite differences approach, nevertheless it can- 
not be applied because the approximations of derivatives 
would not be accurate at discontinuities; even though 
it is common to use finite differences modifying Euler's 
equations with a dissipative term and analyze the limit 
at which such term vanishes 0. Instead, a more ac- 
curate approach used to solve hydrodynamics equations 
considers the use of finite volume methods, which need 
the system of equations to be written in a flux balance 
law form, which in turn requires the definition of conser- 
vative variables as shown below. 

As an illustration of how to write down a balance flux 
equation we construct the first of Euler's equations, the 
one obtained by developing Q for the line element ([!]): 



dtipootau 1 ) H — ^d r (pnaar 2 u r ) = 0, 
dtipoaau*) -I — ^d r (paaar 2 [av r w*]) = 0, 
dt(aD) + \d r {aar 2 v r D) = 0, 



(7) 



which is the first of Euler's equations. The remaining 
equations are obtained in a similar way by developing 
© , which is shown Q . Finally one obtains the following 
set of equations 



d t {aD) + —d r {aar 2 Dv r ) = 0, 

1 2p 

dt(aS r ) H — ^d r (aar 2 [S r v r +p]) = aa — 



'-(S r v r +T+P + D), 



dt(ar) + ^d r (aar 2 (r + p)v r ) = —aa^S r , (8) 
where the set of conservative variables is defined by 



D = Po W, 
S r = pohW 2 v r , 
t = p hW 2 — p — poW. 



(9) 



Then it is possible to write down these equations as a set 
of balance flux type of equations 



d t (au) + -^d r (aar 2 F(u)) = S(u), (10) 

where u is the state vector of conservative variables, F 
the flux vector and S is a source vector given by 





D 




Dv r 


u = 




, F.= 


S r v r + p 




T 




(r + p)v r 



s = 







■2» 



-aa^{S r v r + t + p + D) + aaf 
-aa^Sr 



(11) 



Notice that the term that goes as ~ pjr is singular 
at r = 0. However in order to regularize the equations 
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there, it is possible to split the flux balance form of the 
equations by appropriately splitting the flux vector and 
avoid the presence of such singular term [H Q : 



d t (au) + -^d r {aar 2 ?x{vi)) + 9 r (aaf 2 (u)) = s(u 
where now the fluxes read: 



■). (12) 





Dv r 




" o" 


fl = 


sS r v r 




P 




(t +p)v r 










'■-(S r v r + t + p - 



D) 



(13) 



There is still a usual ingredient when solving problems in 
spherical symmetry, that is, the coordinate singularity at 
r = of the derivative operators in (|f 2j) . This problem is 
solved usually by substituting ^d r f by 3d r 3f for a given 
function /, where now the derivative is with respect to 
r 3 . This result is applied to the second term in equation 



Therefore, we end up with three evolution equations 
for the four variables D, S r ,T,p, or equivalently, for the 
primitive variables po 7 v r ,e^p. This requires to close the 
system, for which an equation of state relating p = p{p) 
is sufficient. In the problems we deal with in this paper 
we choose the gas to obey an ideal gas equation of state 
given by 



P = (r - I )p e, 



(14) 



where T = c p /c v is the ratio between the specific heats, 
sometimes described in terms of the polytropic index n 
such that r = 1 + 1/n [Ioj|. 



B. Spectral decomposition of the spherically 
symmetric relativistic Euler equations 



and their corresponding linearly independent eigenvec- 
tors: 



ri = 



hW(K,~p cj) 
V r 

K 

hW(ti- Pa c1) 



(16) 



1~2 



KW (v r 



hW 
hW 



f r -(A 2 +/3 



-V r (\2+P 

-v r v 



WW) 



(17) 



1 

v r -(\3+P r )/a 



( f rr — v r v r \ 
g^-V(X 3 +^)/a ) 



- 1 



(18) 



where v 2 = g rr v r v r , c s is the speed of sound and (3 r is 
the shift vector, which in our case, according to JT]) is 
zero and g rr = a 2 . On the other hand, a useful property 
of the gas is the speed of sound, which is defined by 



he 2 



X 



\ 



dp 



Op 



(19) 



where for the ideal gas equation of state (|T4")) , \ = e(T — 
1) = p/po an d k = pa(T — 1). Using expression ([3]) 
for the enthalpy, h = 1 + e+ p/po = 1 + f^j^jj 0110 
obtains an expression for the speed of sound in terms of 
the thermodynamical variables 



pr(r - 1) 
p r+p(T-iy 



(20) 



which can be used to calculate the eigenvalues and eigen- 
vectors in (HSJ) and (|16|17|18|) . 



The High Resolution Shock Capturing methods used 
here consider schemes where the spectral elements of the 
Jacobian matrix, A(u) = d , associated to relativistic 
Euler's equations (fTO]) play an important role. Following 
[H, the three eigenvalues of the Jacobian matrix of Euler's 
equations are as follows: 



Ai = -P r + av r 



(15) 



A- 



-/3 r 



1 — v 2 c 2 



v r {l - cl) + v/ C 2(l - v 2 )[g rr (l - v 2 c 2 s ) - v r v r {l - c 2 )} 
A 3 = -P r + 



1 — V 2 C 2 



v r {l - cl) - y/c*(l - v 2 )[g rr (l - v 2 c 2 ) - v r v r (l - C 2)] 



Numerical methods 



Finite Volumes 



Due to the non-linearity of the relativistic Euler equa- 
tions, the presence of discontinuities (shocks) in the hy- 
drodynamical variables is common, even if smooth initial 
data are considered. For this reason, numerical methods 
based in the continuity of the functions, like finite dif- 
ferences, are not suitable to solve this kind of non-linear 
equations. One of the most widely used approaches in the 
literature is the finite volume method; firstly this method 
considers the problem is defined on a mesh of grid points 
that define a cell structure on the space-time, that is, 
the time is restricted to have the discrete set of values 
t n = nAt and the space is considered to be discretized 
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with cells whose spatial centers are defined at Xi = iAx 
as shown in Fig. [TJ Secondly, the numerical methods 
for balance flux type of equations consist in discrctizing 
the equations in their integral form. This is a method 
particularly useful when discontinuities are expected to 
appear as in the present case [f|. 



gfi+1/2 ^ ^ e spatial and temporal average of the sources. 
The key problem consists in the calculation of the tem- 
poral averages of the fluxes which we describe below. 



2. Cell reconstruction 

























AX 










n+1/2 




At 







FIG. 1. In this figure, we present the discretization and cell 
structure of the space-time. Here, the center of cell (j™ +1 / 2 
is located at (t n+1 ^ 2 , xi) and its space-time volume is V — 
AtAx. 

In order to illustrate the finite volume method, con- 
sider the following conservative system of equations in 
one spatial dimension: 



d t u + d x F(u) = S, 



(21) 



where u is the vector of conservative variables, F(u) are 
fluxes which depend on the variables u and S is a vector 
of sources. In this way, the first step to discretize the 
integral form of the equation (pH]) is to take its average 
over a space-time cell c™ +1 / 2 : 



AtAx 
1 



AtA 



Xi+l/2 ft 



,-1/2 
Bi+1/2 



I f x i + l/2 



AtAx 



dtudxdt + 

d x F(u)dxdt 

Sdxdt, 

yn+l/2 



(22) 



where the volume of the cell C^ 1 ' 1 * [ s (t n ,t n+1 ) x 
(%i-i/2> %i+i/2)- Now, as second and final step, by us- 
ing Gauss' theorem, this last equation can be integrated 
to obtain a discretized version of the integral form of the 
system of equations (pH]) : 



u n+1 = u n 



where u™ are the spatial averages of the conservative vari- 

=in+X/2 
i+l/2 



ables, ¥™t}!n are the temporal averages of the fluxes and 



Equation (|23[) provides an evolution rule for the con- 
servative variables, however the difficulty to calculate the 
temporal averaged fluxes F™^ 1 ^ 2 at interfaces between 
cells, still remains. One way to solve this problem defines 
the Godunov-type numerical methods [111 ]. The main 
idea of these methods is to consider a Riemann problem 
at each intercell boundary, which requires to approxi- 
mate the spatial average of the variables u at each cell 
with pieccwisc functions u. 

There are different ways of reconstructing the vari- 
ables. The simplest reconstruction was introduced by 
Godunov and consists in defining the variables to be con- 
stant piecewise @, [ll| . More general reconstructions as- 
sume the variables are linear piecewise, which is the case 
we handle in this paper; in this case, the variables u are 
reconstructed using a minmod slope limiter that restricts 
the slope of the linear functions defining the variables 
within each cell [a, UM '■ 



uf+l/2 = u * + 0i(a»+l/2 - Xi), 

= Q *+l + Vi+l{Xi +l/2 - X i+ x), 



(24) 



where L and R indicate the cell to the left and to the 
right respectively from the intercell boundary we deal 
with. The quantity o\ is calculated as 



minmod(m i _i/2, ^+1/2)- 



(25) 



The function m i+1 / 2 is the derivative of the conservative 
variables u, centered at the cell interfaces 



i+1/2 



Ui+l 



u. 



(26) 



X%-\-\ Xi 

and the minmod slope limiter is defined by 

a if \a\ < \b\ and ab > 
mod(a,b) = { b if \a\ > \b\ and ab > 0. 
if ab<0 

(27) 

The illustration of these reconstructions are shown in 
Fig. [2] The constant piecewise reconstruction provides 
a first order approximation in space of the variables 
whereas the linear reconstruction is second order accu- 
rate. This property impacts on both accuracy and order 
of convergence of the evolution algorithm. In this paper 
of course we use the linear reconstruction which helps 
at achieving convergence of our results. Finally, once 
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we know the functions (|24|) to the right and left of the 
interfaces between cells, we can compute the temporal 
averaged fluxes using an approximate Riemann solver. 




FIG. 2. We show how a function is approximated within each 
cell using Godunov's constant piecewise reconstruction and 
the minmod linear piecewise reconstruction. 



3. Approximate Riemann Solver: HLLE 



4- Calculation of primitive variables 

As mentioned above, the numerical fluxes and sources 
depend both on the conservative and on the primitive 
variables. After each time step within an evolution 
scheme like that in (|23|) during the evolution, one obtains 
new values of the conservative variables D^, S r i, Ti at each 
cell interface i across the grid, then it is required to re- 
construct the primitive variables out of the conservative 
ones in order to account with the necessary information 
to calculate the numerical fluxes and sources (|28|) for the 
expression (fl~3|) . 

In order to do so, first of all we recognize that the vari- 
ables that are being evolved in our system of equations 
are w = ou = (aD,aS r ,ar) := (w\, W2, W3) (see equa- 
tions ©)• From definition ([§]) one can solve for two of 
the primitive variables: 

Pa = w = ! v Vl ~ a2{vr)2 > (30) 

y r = S r = W2 , 3; g 

a 2 (r+p + D) a 2 (wz + ap + wi)' 
whereas the pressure is given by: 



Different approximate Riemann solvers require differ- 
ent characteristic information from the Jacobian matrix, 
for instance the Roe solver (see e.g. pjj]) requires the 
eigenvalues and the eigenvectors. One of the appeal- 
ing properties of the HLLE approximate flux formula is 
that it requires only the eigenvalues of the Jacobian ma- 
trix. The HLLE numerical fluxes formula at the intercell 
boundaries reads 1131: 



A+F(u^ 



-1/2 



+1/2 



- A~F( 



*i+l/2 



+ A+A-(af +1 



P = Poe(r - 1) = (poh - p -p)(T - 1) 
S r D 

P 



= (T-1) 



W 2 V r W 



(r-i; 



D 
W 



S r /v r - DW - pW 2 



+1/2 "f+1/2 



po(r-i) 
po(r-i) 



DW 

t + D(1- W)+p(l - W 2 ) 

DW _ 
w 3 + w 2 {\ -W) + ap(l - W 2 ) 
Wl W 



,(32) 



A+ - A- 



(28) 



where the different As are defined by 



A+ 



1 ■> A 2 ) A 3 > A l 1 A 



max(0, Af , Af , Af 
min(0,Af,A^,Af,Af,A^ 



2 > ^3 . 



)• 



(29) 



Here u L and u R are the values of the conservative vari- 
ables reconstructed at the right and left cells from the 
intercell boundary respectively. Notice however that for 
the cases in this paper, the fluxes and sources in f| 12[) 
depend not only on the conservative variables, but also 
on p and v r which are primitive variables. This requires 
the reconstruction of these primitive variables too. Also, 
the different As depend on the speed of sound, which in 
turn depends on p (see [20]) and is needed at both the left 
and right cells. This is the reason why it is required to 
calculate the conservative and primitive variables. 



where W = W(v r (p)). This defines a trascendent equa- 
tion for p which has to be solved at each cell in the do- 
main. We solve this equation using a Newton-Rapson 
root finder for p at each cell. With this it is possible to 
reconstruct the p in order to obtain p L and p R . Then 
using u L and u R together with (|3"Tj) we calculate v L and 
v R and using (|30|) the rest mass density. Then it is pos- 
sible to calculate the speed of sound on the left and right 
using (f2TJ)) , with this one can calculate the eigenvalues 
of the Jacobian matrix (|15p , which in turn allows one to 
calculate A + and A~ using (|29l) and finally the numerical 
fluxes 



D. Evolution 

The evolution algorithm can be summarized as follows: 

1. Start with given values of the primitive variables 
that contain the physically relevant information of 
the problem. 
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2. Calculate the conservative variables. 

3. Reconstruct all the conservative variables at the left 
and right from the intercell boundaries using (|24|) 
and the pressure by solving (|32j) and then also re- 
constructing to the left and to the right from all 
the intercell boundaries. With this calculate also 
the velocity and rest mass density at the left and 
right cells. 

4. Calculate the speed of sound and the eigenvalues of 
the Jacobian matrix and use them to calculate A + 
and A - . 

5. Use such result to calculate the numerical HLLE 
fluxes IpS]). 

6. Integrate in time the expression (|23p for equations 

o. 

7. Calculate the primitive variables 

8. Repeat from step 3 on. 

In order to evolve the averaged conservative variables 
from time step i" to t n+1 we use the discrete expression 
(|2lfl) with numerical fluxes (|28p corresponding to equa- 
tions (|12[) . which we perform using the method of lines 
(MoL) with a third order TVD Runge-Kutta integrator. 

III. THE BLAST WAVE PROBLEM 

The spherical blast wave problem is a particular real- 
ization of a Riemann problem. Physically, it consists in a 
relativistic gas distributed into two chambers separated 
by a removable spherical membrane located at r = tq. 
Initially the gas in the inner chamber has a higher den- 
sity and pressure than in the outer one, and the velocity 
is zero everywhere. Once the membrane is removed, a 
shock wave moves from the region of higher pressure to 
the region with lower pressure. Also a rarefaction wave 
travels in the opposite direction. Strictly speaking, there 
are various waves propagating in the space, a shock wave 
moving outwards, a rarefaction wave moving inwards and 
between the shock wave and the tail of the rarefaction 
wave two new states are developed which are separated 
by a contact discontinuity. 

We present different physical situations corresponding 
to this problem, based on the numerical solution of the 
relativistic Euler's equations. Specifically, the blast wave 
problem is set on top of the Minkowski space-time. All 
we need to do is set the metric functions to the values 
a = a = 1 in (TTJ). In order to illustrate the physics of 
the spherical blast wave problem, we perform different 
simulations. The parameters we use are shown in Table 
III 

In the weak blast case shown in Fig. the difference 
of pressure is one order of magnitude, whereas in the 
strong blast wave case case in Fig. [H the difference of 



pressure in the initial shock is of three orders of magni- 
tude, respectively. As we can see, the presence of a blast 
wave is more remarkable when the difference of pressure 
is higher, producing velocities close to the speed of light 
and regions where the fluid is supersonic. The case corre- 
sponding to the strong blast wave, has regions where the 
Lorentz factor approaches the value of 4, which indicates 
the relativistic nature of the process. 

An interesting situation is presented in Fig. [S] Due 
to the symmetry of the problem a reverse shock wave 
appears. Unlike the cartesian blast wave case, here the 
two states separated by the contact discontinuity are not 
constant. It happens that in some localized regions the 
pressure is higher to the right side than to the left, pro- 
ducing a shock moving inwards. The velocities reached, 
for this reverse shock wave, are close to the speed of light. 



TABLE I. Initial data configurations for the blast wave prob- 
lem. Subindex "i" is used to represent the initial primitive 
variables of the fluid in the inner chamber of radius 0.5, 
whereas subindex "e" is used to represent the variables in 
the outer chamber r £ [0.5, 1]. 



Case 


Pi 




Pi 




Vi 


V e 


Weak blast 


1.0 


0.1 


1.0 


0.125 


0.0 


0.0 


Strong blast 


133.33 


0.125 


10.0 


1.0 


0.0 


0.0 



IV. THE EVOLUTION OF A TOV STAR 

Concerning the evolution of the TOV star, besides the 
evolution of the fluid as described in detail in the preced- 
ing section, the space-time geometry is allowed to evolve 
according to Einstein's equations, therefore the coupled 
Euler-Eisntein system of equations holds. 

Moreover, in this case the initial data one supplies is 
not arbitrary, instead they have to obey Einstein's field 
equations at initial time. 

In this section we construct various initial configura- 
tions and show how stable and unstable configurations 
behave. 



A. Einstein's equations 

Einstein's equations G Mt , = 87rT^„ for the space-time 
(HJ) and the stress-energy tensor ([2]) in terms of the con- 
servative variables reduce to the following system of 
equations, which coincide with those in [H (see Q for the 
details on the equations): 
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^> 0.3 









0.2 


0.4 


r 


0.6 


0.8 






FIG. 3. Weak blast wave model at t — 0.4 is shown. The initial discontinuity is located at r = 0.5 and the adiabatic index 
is T = 1.4. The spatial resolution used to carry out this numerical simulation is dr — 2 x 10 -4 with a Courant factor of 
At/ Ax = 0.25. Notice that there is small region where the Mach number is bigger than one, which indicates that the fluid is 
supersonic. 



dta = —AirraaS. 
d r a = 
d r a 



Anr(T + D) 



4irr(S r v r + p) + — 



(33) 
(34) 

(35) 



where we have identified the metric function a with the 
mass aspect function using the expression a 2 = 1/(1 — 



2m(r)/r), where m(r) is the mass contained within a 2- 
sphere of radius r. 

Notice that this is an overdctermined constrained evo- 
lution system, that is, the first equation is an evolution 
equation for the metric function a, the second is the 
Hamiltonian constraint and the third equation is a slic- 
ing condition for the lapse a. This system of equations 
allows the evolution of any source provided p, p, D, S r 
and r, however in order to represent a solution of Ein- 
stein's equations they need to satisfy such equations at 
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60 
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FIG. 4. Strong blast wave model. The numerical parameters are as in the previous case and the gas obeys the same equation of 
state with F = 1.4. In this case the Lorentz factor reaches values near 4. The region where the velocity of the fluid is supersonic 
is wider and the velocity there higher than in the previous case. 



initial time, which wc describe next. 



B. The initial value problem 

A TOV star is described as a spherically symmetric, 
static system that obeys Einstein's equations sourced by 
a perfect fluid that obeys a polytropic equation of state. 

In order to solve the initial value problem we assume 
the space-time metric is static and momentarily will use 
m(r) instead of a in order to maintain the standard no- 



tation for the construction of TOV stars (see e.g. 
Then we start with the line element ([IJ rewritten as: 



ds 1 



-a(r) 2 dt 2 



dr 2 



1 



2m(r) 



'-de 2 



: sin 2 9d(j) 2 , (36) 



where we have assumed the system is time-symmetric at 
t = and the metric functions and the gas functions 
depend only on r. We also assume the gas obeys initially 
a polytropic equation of state p = Kp T . Using [l4[ one 
arrives at the following set of equations: 
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FIG. 5. In this figure, the presence of a reverse shock wave, in the spherical blast wave problem, is shown. The parameters we 
use, in this example, are the following: domain of the simulation r £ [0, 12], numerical resolution dr = 0.004, Courant Factor 
0.25, adiabatic index 1.4, initial discotinuity ro = 3.0 and the initial primitive variables are pi — 13.33, p e = 0.1, pi = 10, p e = 1 
and Vi — v e = 0. 
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dm 2 
— — = 47rr p, 
dr 



dp 
dr 



1 da 
a dr 



r(p + p) 





47rr 3 p\ 


( 2m\ 


I 1 --) 




m ) 





-(p+p) 
1 



4irr 3 p 



r(r — 2m) ' 
dp m + Airr 3 p 



p + p dr r(r — 2m) 



(37) 

(38) 
(39) 



where we have used p = po(l+ e ) an d Poh = po(l+e)+p = 
p + p. This system of ordinary equations constitutes the 
conditions a TOV star satisfies at initial time, and has to 
be integrated outwards from r = up to r = r max . We 
solve these equations using a fourth order Runge-Kutta 
integrator on top of the same grid defined for the fluid 
equations described above. The initial conditions for the 
integration of the variables are: 

1. m(0) = 0, because the integrated mass up to 
there is zero. Another interpretation is that the 
gravitational field at the origin is zero, and thus 
a(r) = 1 corresponds to the flat space, which im- 
plies m(0) = 0. 

2. p(0) = Kp T cl where p c is the central value of the 
rest mass density. In the whole domain it happens 
that on the one hand po = [p/K) l / v and on the 
other, from an ideal gas equation of state p = {T — 
l)poe, =S> p e = p/(T-l); therefore p = po{l+e) = 
Pa + Pot = (p/K) 1 ' F + p/(r — 1) is the source of 

(EH). 



3. a(0) = ao is an arbitrary given initial central value 
for the lapse. Notice in (|3"9")l that the solution 
can be rescaled multiplying by a constant, which 
preferably will be chosen such that at the numeri- 
cal boundary satisfies a(r max ) = l/a(r max ), which 
is a condition that Schwarzschild's solution satisfies 
and we expect to happen at r = r max . 

4. The value of p c turns to be the input parameter 
that determines the configuration, and corresponds 
to the central value of the total energy density. The 
result is that for each value of p c a configuration can 
be constructed. 

Two observations are in turn. The first one concerns 
the point r = 0, because there equations (|3"5)) and (|3"9"|) are 
singular. What is usually done is to Taylor expand the 
singular factor and get approximate equations for small 
values of r: 

m + 47rr 3 p 
r 2 — 2mr 

m(0) + m'(Q)r + |m"(0)r 2 + im"'(0)r 3 + 0(r 4 ) + 4?rr 3 p 
~ r 2 - 2r(m(0) + m'(0)r + ±m"(0)r 2 + ±m"'(0)r 3 + 0(r 4 )) 
_ Awpr/3 + Anrp 

~ 1 ~8npr 2 /3 [ ' 



where equation (|3T|) was used to calculate the deriva- 
tives of m: dm/dr\ r= o — Ai:r 2 p\ r=Q = 0, d 2 m/dr 2 \ r=0 = 
4n(2rp + r 2 p)\ r=Q = and d 3 m/dr 3 \ r=0 = 4ir(2p + 2rp + 
r 2 p)| r= o = 8np. Then equations (|38|) and (|39|) are ap- 
proximated for small r by 



dp . 4ir pr / 3 + 4irrp 

= -(P + P) l _ 8lTpr 2/ 3 ' 

1 da Anpr /3 + Airrp 
a dr 1 — 8irpr 2 /3 



(41) 



These approximate regular equations are the ones to be 
programmed for at least the first mesh point located at 
r = Ar where Ar is the spatial resolution of the mesh. 

The second observation is related to the divergence of 
the specific enthalpy ([3]) when po approaches zero, which 
in theory would happen from the star's surface to infinity 
where there is only vacuum. It is usually set an external 
atmosphere, that is, a minimum value is assumed for po 
than can be hidden within numerical errors and allows 
the convergence of the numerical calculations, however it 
happens to be a mere numerical artifact at the moment 
and as far as we can tell, there is no theory behind the ap- 
propriate value of the atmosphere density p a tm = floor. 
The value of floor rather depends on the specific problem 
to be solved. 

Considering this ingredient one can define the radius 
R of the TOV star as the minimum radius r = R at 
which po = floor. On the other hand, the total mass 
of the TOV star is Mr = m(R), whereas the rest mass 
of the star is the spatial integral of po given by Mo = 
47r J Q R por 2 a(r)dr. The difference between the total and 
the rest mass of the star determines whether or not the 
system is gravitationally bounded. 

As an example of a TOV star configuration we show 
in Fig. [S]the functions for T = 2, a polytropic constant 
K = 1 and a central density p c = 0.42. 

The result of integrating the TOV equations for vari- 
ous values of p c is summarized in Fig. [71 where we plot 
the total and rest mass for two different classes of equa- 
tions of state, r = 2 and T = 5/3 for several values of 
p c . Each point in the curves corresponds to a value of p c 
and therefore defines a TOV star configuration. The first 
plot corresponds to an ultrarelativistic case whereas the 
second serves to model a fcrmionic gas and is a simple 
approximate model of white dwarfs. The maximum in 
the plots indicates the threshold between stable and un- 
stable configurations, that is, configurations to the left 
of the maximum oscillate under perturbations whereas 
those to the right collapse and form black holes if they 
are perturbed, because these systems are gravitationally 
bound since Mq > Mt for the values of p c shown. We 
also indicate in Fig. [7] four particular configurations, two 
stable and two unstable that we use to illustrate their 
evolution. 
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K=1, r=2, p c =0.42 




K=100, T=2 



FIG. 6. Initial data as found for T = 2, (n = 1), p c = 0.42, 
a c = 0.5, If = 1. The result is as follows, the total mass 
M = 0.1616, the rest mass M = 0.177, R = 0.7045. The 
lapse a has been rescaled such that a(r ma x) = l/a(r mtJ x) as 
expected to happen for Schwarzschild's solution. 



Summarizing, the information required to start up the 
evolution of a TOV star has now been calculated at each 
cell, and is the following: 



a 
a 
P 

Po 

W : 



h 
D 

T 



(numerically integrated), 
(numerically integrated), 
(numerically integrated), 

0, 

1 



P 

Po(r-i)' 

1 + e + p/p , 
PoW, 

Po hW 2 a 2 v r , 
p hW 2 - p- p W. 



Then the evolution of the system is ruled by the Einstein- 
Euler system as described next. 



C. The evolution 

The system of equations is the one composed of Euler's 
equations (|12[) and the overdetermined system of Ein- 
stein's equations (|33H35j) . The whole system is started 
with the initial data corresponding to a TOV star. 
Among Einstein's equations we choose to solve (f3"3"f for a 



and the remaining (|34j) is the Hamiltonian constraint we 
use to monitor the evolution. Notice that the equation 
for a is an ODE in r, that we integrate every time step 
during the evolution. 




0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 



POC 

K=10, T=5/3 




Total mass 

Rest mass 

Critical point • 
Special cases ■ 



0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 

POc 



FIG. 7. Mass vs central density diagram for the cases 
K = 100, r = 2 and K = 10, T = 5/3. The filled circle 
indicates the location of the maximum possible mass and also 
the threshold between the stable and unstable brances. Those 
configurations to the left of the maximum are stable and those 
to the right are unstable. The fact that the rest mass has a 
bigger value than the total mass indicates that the system is 
gravitationally bounded. This eventually implies that those 
configurations belonging to the unstable branch should col- 
lapse and form black holes. Configurations marked with a 
filled square correspond to particular configurations we evolve 
to illustrate the different behaviors of stable and unstable con- 
figurations. 



Thus the algorithm for the evolution provided given 
values of K and T is as follows: 

• Construct a TOV configuration as initial data for 
a given central density p c . 

• Use the evolution equations and calculate new D, 
S r , t and simultaneously integrate in time the mo- 
mentum constraint (|33p for a at every time step. 

For each intermediate time-step of the MoL 
integrator 

— Express the primitive variables in terms of the 
conservative variables and reconstruct to the 
left and right from intercell boundaries the 
values for p, v r , po, and the conservative vari- 
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K=100, r=2, p 0c =0.001 

1.1722 | 1 1 1 1 
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4.4 




f 



K=10,r=5l3, p 0c =0.0006 

1.1284 | 1 1 1 , 
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FIG. 8. Maximum of a(r) for the stable cases K = 100, r = 2 and p c = 0.001 and K = 10, T = 5/3 and p c = 0.0006. The 
metric function responds to the perturbation due to numerical errors and its maximum oscillates. The constant line indicates 
the value calculated at initial time for the maximum of a which should be maintained constant in the continuum limit. Also 
shown is the convergence factor of the L\ norm of the Hamiltonian constraint of Einstein's equations. This simulations use 
6000 cells in a domain r £ [0, 500] for the V = 2case and 3000 cells for the case F — 5/3 in a domain r 6 [0, 500]. In both cases 
the atmosphere rest mass density is floor = 10 -13 . 



ablcs D, S r , t in order to obtain the necessary 
information to construct the numerical fluxes. 

— Apply boundary outflux conditions to the con- 
servative variables and extrapolate for a. In 
our conservative formulation it requires only 
to copy the values of the conservative variables 
at boundaries from the point next to it. 

— Integrate ([33]) for a, and rescale it such 
that at the boundary it satisfies a(r max ) = 
l/a(r max ). 

• After a full time step, calculate the Hamiltonian 
constraint in order to monitor the convergence 
of the results. 

Let us explain what should happen during the evo- 
lution. In the continuum limit the TOV configurations 
should remain time independent all the way, because they 
are static solutions to Einstein's equations. If a pertur- 
bation is applied (for example, a small amplitude shell 
pulse added to the density), the geometry and matter 
quantities would oscillate around the equilibrium values, 



whereas an unstable configuration would collapse and 
form a black hole. Nevertheless, we are using numerical 
methods and as shown above, all our calculations involve 
an intrinsic error. We then take advantage of such fact 
and use such error as the perturbation of the equilibrium 
configurations. Therefore stable configurations would os- 
cillate around the equilibrium values, whereas unstable 
configurations eventually will collapse due to a perturba- 
tion triggered by the numerical errors. 

In order to illustrate the evolution of TOV stars we 
choose two stable configurations and show some results 
in Fig. [5] On the one hand we show the maximum of 
the metric function a in time which shows a periodic os- 
cillation. The reason is that we are solving numerically 
the initial value problem, and also we are integrating nu- 
merically with a finite accuracy, then there is a numerical 
error introduced in our calculations at initial time which 
works as a perturbation whose effects converge to zero in 
the continuum limit fl5| . What is more important is that 
the metric function remains nearly time-independent, as 
expected for a stable equilibrium configuration. On the 
other hand we show the convergence of the Hamiltonian 
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K=100, T=2, p c =0.004 K=100, r=2, p c =0.004 
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r r 

FIG. 9. We show snapshots of the metric functions for the unstable cases indicated in Fig. [7] (K — 100, T = 2, p c = 0.004) 
with total mass Mt = 1.623 and (K = 10, T — 5/3, p c = 0.0025) with Mt = 1.475. The lapse collapses to zero with time, 
which indicates that an apparent horizon has formed, and in turn implies that external to such apparent horizon there is an 
event horizon; observe that the lapse approaches zero until r ~ 2Mt- Notice also that the metric function a diverges at a 
similar location of the horizon radius, which is an effect of the slice stretching that occurs during a black hole formation in 
non-penetrating coordinates. 



constraint, which is necessary to verify that we are truly 
solving the full set of Einstein's equations (remember 
that the Hamiltonian constraint (|34p is not being solved, 
only monitored). We are verifying the convergence of 
our results by doubling the resolution, which means that 
a convergence factor is defined as 2^ where Q is the or- 
der of convergence |15j . Then from Fig. [5] we know our 
results converge within order 1.6 and order 2, which is 
consistent with the approximations we have made in all 
the methods used. 

We also show the evolution of two unstable configura- 
tions in Fig. 121 In this case the metric functions a and 
a do not remain nearly time independent as in the sta- 
ble cases, where snapshots of the metric functions would 
be seen as a single curve. Instead, the lapse collapses to 
zero in a localized region, which in the coordinates we use 
means that an apparent horizon and that a black hole has 
been formed. Also the function a diverges near the loca- 
tion of the horizon, which is due to the slice stretching 
effect of the normal coordinates we are using. 

In order to make sure that a black hole has formed 
we track a bundle of outgoing null geodesies starting at 



about t ~ 103, which we show in Fig. [TU] for one of the 
collapsing configurations. The null rays shown indicate 
the behavior of null 2-sphcrcs. Near the event horizon 
these null surfaces should diverge toward the singular- 
ity and toward future null infinity. We show the null 
geodesies until our simulation remains accurate, which 
happens until the aforementioned problem of the coordi- 
nates we are using appears. However this small window 
in time allows one to appreciate the divergence of the null 
spheres and thus infer that the event horizon is contained 
into the set of null rays shown. This guarantees that not 
only an apparent horizon has been formed through the 
gauge dependent condition a ~ but also an event hori- 
zon, which is gauge independent. It is possible to see that 
the event horizon grows due to the accretion of the gas 
and tends to stabilize at a radius nearly twice the mass 
of the initial configuration. 
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Event horizon of a collapsing configuration 
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FIG. 10. We show a bundle of null rays, which due to the 
symmetry and coordinates we are using, represent the behav- 
ior of 2-dimensional null spheres. What we show is that some 
of these rays diverge toward the singularity and other will 
diverge outwards toward infinity, and more precisely toward 
future null infinity [13]. The growing null sphere in between 
the two behaviors would be the event horizon. In theory, an 
event horizon would be well defined if we guarantee that out- 
going null rays diverging outwards reach future null infinity, 
however we only show a small chunk of the space-time where 
we can measure the divergence of these null surfaces due to 
the aforementioned slice stretching drawback of the coordi- 
nates we are using. 

V. FINAL COMMENTS 

We have shown in detail a particular sort of implemen- 
tation of numerical rclativistic hydrodynamics solutions, 
of spherically symmetric cases in spherical coordinates. 
The steps specified in the paper are also useful for differ- 
ent choices of numerical approximations described here. 

Specifically, related to the treatment of relativistic hy- 
drodynamics, we only use a particular flux formula for 
the numerical solution of the Riemann problems at the 
intercell boundaries. There are several other choices like 
the Roc, Marquina, HLL, HLLC, etcetera flux formulas. 



Also, for the cell reconstruction of variables, other choices 
aside the minmod limiter are well studied like the MC 
(linear monotonic centered), PPM (parabolic piecewise 
method), etc. 

For instance in [§[ the authors use the combination 
HLLE flux formula and MC limiter for the evolution of 
TOV stars, and in [n| the evolution of TOV stars is 
implemented using Marquina and Roe fluxes with MC 
and minmod slope reconstructors. 

We found appropriate to choose a single combination of 
numerical methods in order to be as specific and detailed 
as possible. 

We also want to mention other aspects inherent to 
these numerical methods. Particularly interesting is that 
the enthalpy diverges when the rest mass density ap- 
proaches zero ([3]), and the implementation of an atmo- 
sphere is necessary. However, so far there is no theory or 
explanation about what values of the atmosphere den- 
sity are to be applied, and in the best cases (as here) 
convergence tests are used to support the numerical re- 
sults, and the values used for such external density is 
justified as long as the numerical results in terms of ac- 
curacy and convergence are achieved. A potential recipe 
for an atmosphere with a different equation of state may 
ameliorate this problem Q. 

Even though the density at the atmosphere is small, 
the fluid may develop highly rclativistic speeds, which 
eventually may produce intractable shocks at the star 
surface. Therefore the atmosphere requires a rather ad 
hoc treatment, like artificial limitation to the speed of the 
fluid at the atmosphere, etc. In particular in our simula- 
tions of TOV stars we have used this type of condition. 
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